subroutine rks1
use var
implicit none

    integer i,m,j
    real(8) du,dv
    real(8) Qmean_n_k(1:gridsx,1:gridsy,5)
    
    
    
    call comput_flux1
    call comput_flux2    
    call rhs(Qmean,Qmean_n_k)
    do j=1,gridsy
        do i=1,gridsx
            Qmean(i,j,1)=Qmean_n(i,j,1)-dt*Qmean_n_k(i,j,1)
            
            Qmean(i,j,2)=Qmean_n(i,j,2)-dt*Qmean_n_k(i,j,2)
            
            Qmean(i,j,3)=Qmean_n(i,j,3)-dt*Qmean_n_k(i,j,3)
            
            Qmean(i,j,4)=Qmean_n(i,j,4)-dt*Qmean_n_k(i,j,4)
            
            !write(*,'(I4,2X,100(ES12.4,1X))')i,j,Qmean(i,j,4)
            !do m=1,4
            !    du=-(flux1(i,j,m)-flux1(i-1,j,m))/dx
            !    dv=-(flux2(i,j,m)-flux2(i,j-1,m))/dy
            !    Qmean(i,j,m)=Qmean_n(i,j,m)+dt*du+dt*dv
            !end do            
            !Qmean(i,j,5)=Qmean_n(i,j,5)-dt*Qmean_n_k(i,j,5)        
        enddo         
    end do
    call comput_original_var
    call boundary
    
    
    
    
    

end subroutine 